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Abstract 


A novel  analysis  of  the  steady-state  probabilities  of  a class  of  infinite 
Markov  chains  is  given.  Markov  chains  of  this  type  appear  in  the  study  of 
bulk  queues  and  a variety  of  other  stochastic  models.  Algorithms,  which 
involve  only  real  arithmetic  and  avoid  the  traditional  analysis,  based  on 
Rouche's  theorem,  are  presented. 
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probability,  Bailey's  bulk  queue,  Moran  dam 
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1.  Introduction 


There  is  a long  tradition  in  queueing  theory  of  presenting  the  probability 
distributions  of  interest  in  the  form  of  Laplace-Stielt jes,  Fourier  or  generat- 
ing function  transforms.  These  transforms  often  involve  several  unknown 
parameters,  which  need  to  be  evaluated  by  the  solution  of  an  auxiliary  system 
of  linear  equations  with  complex  coefficients.  The  coefficients  of  that 
auxiliary  system  depend,  in  addition,  on  the  roots  of  a (usually)  transcendental 
equation  inside  some  region  of  the  complex  plane.  The  numerical  computation 
of  the  probability  distributions  of  interest,  by  this  classical  method  based 
on  an  application  of  Rouche's  theorem,  is  in  practice  quite  difficult.  It  can 
also  be  numerically  hazardous,  when  the  roots  of  the  transcendental  equation 
lie  close  together  or  coincide.  From  an  aesthetic  viewpoint  it  is  unattractive, 
as  it  involves  several  steps  without  a clear  probabilistic  significance.  In 
spite  of  its  use  in  literally  hundreds  of  papers  on  queues,  this  method  has 
also  been  largely  ignored  by  the  practitioner  who  attacks  the  numerics  of 
queueing  problems  by  ad  hoc  truncation  procedures  or  by  computer  simulation. 

It  is  possible  however  to  show  that  many  queueing  problems  can  be  solved 
by  a purely  probabilistic  method  which  requires  little  or  no  complex  analysis. 
This  approach  leads,  at  least  in  the  range  of  practical  interest,  to  efficient 
and  stable  algorithms  involving  only  real  arithmetic.  The  theoretical  back- 
ground material  for  this  procedure  is  discussed  in  Neuts  [8,  11]  and  involves 
some  elementary  functional  analysis.  Applications  to  the  queue  with  semi- 
Markov  service  times  and  to  an  M/G/l  queue  with  several  types  of  customers  are 
available  in  Neuts  [12,  13].  Extensive  numerical  examples  and  computer  programs 
are  appended  to  the  preprint  of  the  paper  [13]. 

In  the  present  paper,  which  is  partly  expository,  we  shall  begin  by  giving 
a detailed  discussion  of  the  class  of  Markov  chains  with  a transition  probabil- 
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ity  matrix  P of  the  form  given  in  Formula  (1).  Chains  of  this  type  arise 
frequently  in  queueing  models.  In  the  latter  part  of  the  paper,  we  shall  show 
how  a large  variety  of  explicit  or  computationally  tractable  new  results  can 
be  obtained  for  some  specific  models,  discussed  in  the  literature. 

We  consider  Markov  chains  with  transition  probability  matrix  P,  given  by 
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respectively.  Their 

probability  generating  functions,  used  only  for  analytic  convenience,  will  be 
denoted  by  A*(z)  and  B*(z),  0 < k < m - 1,  respectively  and  in  order  to  avoid 
triviality,  we  assume  that  m > 1. 


A.  Trreducibility  and  Aperiodicity 

In  applications,  it  is  usually  obvious  from  the  specific  form  of  the 
entries  that  the  chain  P is  irreducible  and  aperiodic,  but  a general  discussion 
of  these  properties  is  tedious  and  not  very  interesting.  It  is  easy  to  see 
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that  under  the  two  conditions,  stated  below,  the  chain  is  both  irreducible  and 
aperiodic. 

Condition  1.  Let  d = g.c.d  {k:  a^  > 0,  k > 1},  then  we  assume  that 
(2)  g.c.d. (m,d)  = 1,  and  a^  > 0. 


Condition  2.  Let  the  matrix  be  defined  by 


(5)  B. 
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Much  of  the  general  theory  remains  valid  for  stochastic  matrices,  which 
can  be  written  in  the  partitioned  form  (4),  but  in  the  present  case  the 
special  form  of  the  matrices  , j >0,  leads  to  some  simplifying  properties, 
noted  in  Lemma  1 below. 

We  first  introduce  some  notation.  The  stochastic  matrices  A and  B are 
defined  by 


(6) 


A = 


Furthermore,  with  e.  = 
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and  we  denote  the  invariant  probability  vector  of  A by  _rr. 
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Lemma  1 


The  matrix  A is  a circulant  and 


— m — 


If  U is  an  integer-valued  random  variable  with  density  {a^},  then 


0 W-k-1, 

a.  = E[ — - — ], 
k m 


for  1 < k < m, 


where  [•]  is  the  integer  part  function.  Furthermore  the  inner  product  tt  £ 


is  given  by 


0 * 
it  a = m a*. 


Proof. 


It  is  clear  that  A is  a circulant  and  this  implies  (8).  See  e.g. 
Rosenblatt  [16]  p.  52.  Furthermore 


0 ~ „rU+k-l, 

a.  = 2,  v a , i = E[ ], 

k v=l  r=0  mv-k+r+1  m 


for  1 < k < m.  For  every  integer  N > 0,  we  have  the  elementary  equality 


U1  X1  , . 1 

•ll  = N, 

J-i  m 


which  clearly  implies  (10). 


Henceforth  we  shall  refer  to  the  set  of  states  {mr,  mr+1,  ...,  mr+m-1) 
as  level  r,  for  r > 0.  For  convenience  we  shall  refer  to  the  original  state 
mr+j-1  also  as  state  (r,j),  for  r > 0,  1 < j < m. 


2.  A First  Passage  Problem 

In  [8],  we  examined  the  first  passage  from  the  state  (r+1,  j)  to  the 
level  r.  By  G^,(k),  k > 1,  we  denoted  the  probability  that,  starting  in  the 
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state  (r+l,j),  the  level  r is  reached  for  the  first  time  after  exactly  k 
transitions  by  a visit  to  the  state  (r,j').  The  sequence  of  matrices  G(k)  = 
tGj j , (k) } , k > 1,  is  of  basic  importance.  We  showed  that  its  matrix  generating 
function 


G(z)  = kl±  G (k) z , 


satisfies  the  matrix  functional  equation 


G(z)  = E zA  G (z), 
v=0  v ’ 


for  z < 1. 


In  an  appropriately  defined  set  of  matrices  G(z)  is  also  the  unique  solution 


to  Equation  (14)  for  0 < z < 1.  The  following  theorems  summarize  results 
proved  in  [8]  and  [11]. 


Theorem  1 

Defining  the  corresponding  first  passage  probability  matrices  G^(k), 
k > i,  between  level  r+i  and  level  r,  the  matrix  generating  function  of  the 
sequence  {G^(k)},  is  given  by  the  i-th  power  G*(z)  of  G(z)  , for  i > 1. 

The  matrix  G,  defined  by 


G = | G(k)  = lim  G(z) , 

z-*-l- 


is  stochastic  if  and  only  if 


0 -1  * 

p=ira  = ra  a*  < 1. 


The  matrix  G is  then  the  unique  solution  of  the  equation 


oo  V 

G « LAG, 
v=0  v 


in  the  set  of  substochastic  matrices.  The  minimal  solution  G to  (17)  in  the 
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set  of  substochastic  matrices  is  strictly  positive.  If  p > 1,  G is  strictly 
substochastic  and  for  p < 1,  we  have  that  G = G is  stochastic. 

Remarks 

a.  The  numerical  solution  to  (17)  was  examined  in  111].  For  most 
practical  purposes,  G may  be  calculated  efficiently  by  successive  substitutions 
in  the  modified  form 

(18)  C - (I  - A^'1  JQ 

V*1 

of  Equation  (17),  starting  with  G = 0. 

b.  The  special  structure  of  the  matrices  A^,  v > 0,  does  not  appear  to 
lead  to  any  tractable  special  properties  of  the  matrix  G.  In  particular  G is 
not  a circulant.  This  is  shown  by  the  example,  where  m = 2 and 


and  A^  = 0,  for  v > 3.  This  rare  tractable  example  has  the  solution  matrix  G, 
given  by 

|(/T7  - 3)  |(5  - vT7) 

2/l7  - 8 9 - 2vT7 

For  p < 1,  the  stochastic  matrix  G has  an  invariant  probability  vector  £, 
which  satisfies  £ = £ G,  and  £<5  = 1.  By  G we  denote  the  m x m matrix  with 
m identical  rows  given  by  £.  It  is  well-known  that  G and  G commute  and  that 
the  matrix  I - G + G is  nonsingular. 

We  define  the  matrix  M by 

M = G’ (1-)  = E vG(v) , 


(19) 
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and  the  vector  y by  p = M e_.  The  matrix  AOj  ) is  the  diagonal  matrix 
diag(aG,  a^) . The  following  useful  results  were  established  in  [11]. 


Theorem  2 

When  p = 1,  the  matrix  M is  infinite.  When  p < 1,  the  matrix  M is  the 
solution  of  the  equation 


(20) 


oo  1 v i-v-1 

M = G+  .L  A.  In  G M G 
1=1  i v=0 


The  vector  p satisfies 


(21) 


i-1  v 

£ + ill  Ai  v=0  G 


and  is  given  explicitly  by 

(22)  p = (I-G+G) [I-A+G-A(a°)G]_1e. 


The  scalar  product  £ _p  is  given  by 
(23)  £ p.  = (l-p)  " = m(m-a*)  1. 


Remarks 

a.  In  the  sequel  we  shall  express  a large  number  of  features  of  the 
Markov  chain  P in  terms  of  the  matrix  G and  the  vectors  £ and  p_.  These 
entities  are  the  crucial  elements  to  be  evaluated  in  a numerical  solution  of 
these  features.  The  equation  (23)  is  a powerful  check  on  the  numerical 
accuracy  of  the  computed  values  of  and  p_. 

b.  The  entries  of  the  matrices  G and  M and  of  various  other  matrices 
expressible  in  terms  of  these  have  a number  of  interesting  interpretations 

in  the  queueing  problems,  which  lead  to  Markov  chains  of  the  type  P.  Some  of 


these  interpretations  will  be  discussed  below. 
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3.  The  Invariant  Probability  Vector  of  the  Matrix  P 

The  invariant  probability  vector  21  of  the  matrix  P satisfies  the 
stationarity  equations 


(24) 


x . 
3 


m-1  j 

x b . + , ln  x , . . a.  , 
v=0  v vj  k=0  m+j-k  k 


for  j > 0. 


The  probability  generating  function 


(25)  X(z)  = kf0  xkzk, 


for  |z(  < 1, 


is  given  by 


(26) 


X(z)  = 


% xvzV[AMz)  - 
A*(z) 


z B*(z) 


We  see  that  the  right-hand  side  depends 


jj_  = (Xq,x^,  . . . .x^  ^)  , which  we  shall  now  determine  by  a probabilistic  argument. 
The  matrix  K is  defined  by 


(27) 


K = Zn  B G , 
v=0  v 


and  is  clearly  stochastic. 


Lemma  2 


The  matrix  K is  irreducible. 


Proof 

If  K were  reducible,  its  rows  and  columns  could  be  permuted  in  the  same 
way,  so  that  K could  be  written  in  the  form 


' K, 


K„ 


(28) 


K = 


10 


This  implies  that  the  matrices  and  B^,  for  v > 1,  must  be  of  the  form 


[29) 


B0  =l 


Bo(1)  Bo(2) 


V3)< 


BV(D  Bv(2) 


, for  v > 1, 


since  the  matrix  G is  strictly  positive. 

This  further  shows  that  Bq(3)  is  stochastic,  and  therefore  that  I - B( 
is  singular.  This  contradicts  Condition  2,  which  we  assumed  to  hold. 

The  invariant  probability  vector  of  K is  therefore  well-defined  and 
satisfies 


:30)  I = IK,  1 e = 1. 

We  further  define  the  vector  by 

oo  V—  1 -i 

31)  0 = e + I,  B .ln  GJy, 

— — v=l  v j=0  — 

where  G and  y are  as  defined  in  Section  2. 


Lemma  3 

The  vector  ! is  given  explicitly  by 
(32)  0_  = e + (1-p )-16°  + (B-K)  (I-G+G)_1jj. 

Proof 

Since  GV(I-G+G)  = GV  - G^1  + G,  for  v > 0,  we  obtain 

;33)  v|L  B^  Gj(I-G+G)  = B - K + A(!°)G, 

where  A(B°)  = diag(6°, . . . ,6°) . Since  I - G + G is  nonsingular  we  obtain 
— 1 m 

9.=  e+  [B-K+A(_6°)G ] (I-G+G)_1y 

= e + (1-p)"1!0  + (B-KHI-G+G)"1!, 


34) 
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by  noting  that  A(3°)G  £ = (g.  £)_B°  = (1-p) 

Corollary  1 

The  scalar  product  y_  £ is  given  by 

(35)  2 6 = 1 + (l-p)"1!  6°  - x(I-B)(I-GfG)-1jj. 

Proof 

From  Formula  (32)  and  the  definition  of  jr- 
Theorem  3 

The  vector  £ = (xn,...,x  .)  is  given  by 

U m-i 

(36)  1=  (Xi)"1!- 
Proof 

We  consider  the  successive  visits  to  the  set  of  states  (0, 1 , . . . ,m-l } 

in  the  Markov  chain  P.  If  J denotes  the  state  at  the  n-th  visit  and  X 

n n 

the  time  between  the  (n-l)st  and  the  n-th  visit,  then  the  sequence  of  pairs 

{(J  ,X  ),  n > 0}  with  X„  = 0,  is  easily  seen  to  be  a Markov  renewal  sequence 
n n - 0 

with  m states.  The  transition  probability  matrix  is  of  lattice  type.  An 
easy  first  passage  argument  shows  that  the  probability  generating  matrix 
function  K(z)  of  that  transition  probability  matrix  is  given  by 

(37)  K(z)  - jf0  z BvGV(z),  for  |z|  < 1. 

We  clearly  have  that  the  matrix  K,  defined  in  (27),  and  the  vector  £, 
defined  in  (31),  satisfy 

(38)  K.  = K(l-),  £ = K’ (l-)e. 


Application  of  a classical  theorem  on  Markov  renewal  processes,  see  e.g. 


I 


^inlar  [2],  Thm.  6.12,  p.  155,  or  Hunter  [4],  Thm.  2.11,  p.  196,  yields  that 
the  mean  recurrence  time  of  the  state  (0,j)  in  the  finite  Markov  renewal 
process  is  given  by 


Ej  ■ w).;1, 


for  1 < j < m. 


Since  is  also  the  expected  number  of  transitions  between  successive 
returns  to  the  states  j-1  in  the  Markov  chain  P,  we  have  that 


V<Vi)  • 


for  1 < j < m. 


but  this  is  equivalent  to  (36). 


Corollary  2 


m_ 1 / n \ 1 

vi0  % = (I  0)  . 


and  the  right-hand  side  is  explicitly  given  by  applying  Formula  (35). 

The  higher  terms  x^. , j > m,  of  the  stationary  probability  density  may  be 

computed  recursively  once  the  initial  terms  x_,  x, , ...,  x , are  known. 

0 1 m-1 

From  the  stationarity  equations  (24)  , we  obtain 

(42)  x = a ^[x.  - xb  -.1.  x , ,.  a ],  for  j > 0. 

m+j  0 j v=0  v vj  k=l  m-k+j  k - 

On  computers  with  a short  word  length,  the  recurrence  relation  (42)  is  for 
large  values  of  j quite  sensitive  to  loss  of  significance,  because  there  are 
a substantial  number  of  small  terms  to  be  subtracted  from  the  small  positive 
value  of  Xj . For  the  case  of  the  M/G/l  queue.  Dr.  Paul  J.  Burke  suggested  to 
the  author  an  alternate  form  [10]  of  the  recursion  formula  for  the  x ^ , which 
is  far  less  sensitive  to  loss  of  signifcance.  For  the  Markov  chain  P,  this 


alternate  form  is  given  by  the  following  result. 
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Theorem  4 


The  quantities  x^ , for  j > m,  may  be  recursively  computed  by  use  of  the 


formulas 


(43) 


m-1  - m-1 

x * a_  [ x b - £.  x ] , 

m 0 v=0  v vO  v=0  v 

-l  j - m-1 

x . . - an  { £,  x . . a + £~  x b . 
nrt\j  0 v=l  nrt-j-v  v v=0  v vj 


n-1 


j-1 


- [ £.,  x + £n  x , ] } , for  1 < j < m-2, 

v=j+l  v v=0  m+v  * - J - 

_i  j « m-1  - m-1 

x , . = an  [ £,  x . a + £_  x b . - £ x ] , for  j > nrt-1, 

m+j  0 v=l  m+j-v  v v=0  v vj  v=l  m+j-v  - 


where 


(44) 


bvj  = 1 " k=0  bvk’  for  0 < v < m-1,  j > 0, 


aj  = 1 ' kio  V 


for  j > 0. 


Proof 


The  generating  function  X^(z)  = z m[X(z)  - x^z'],  satisfies 

m m—  1 \»  m 1 

(45)  z\(z)  + v|0  xvz  = v|0  xyB*(z)  + A*(z)X1(z), 


m-1 


and  this  leads  to 


(46) 


1-z 
1-z  "1 


•X,(z)  = 


1-A*(z)  , . . m;1  1 Bv(z)  _ m:1  l-zV 

1-Z  1 Z v=0  XV  1-z  v=l  xv  1-z 


Recalling  that 
(47) 


l~A*(z)  _ “ 2 ,k 

1-z  " k=0  V 5 


1-B*(z)  . . 

v _ » ; k 

1-z  k=0  vk 


for  0 < v < m-1,  and  expanding  both  sides  in  (46),  we  obtain  upon  equating  the 
coefficients  of  zJ , the  expressions  given  in  Formula  (43). 
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Corollary  3 


The  equality 


(m-a*) [ 1 - Ex] 
v=Q  v 


vi0  Xv(B*'v)* 


holds. 


Proof 


By  letting  z -*  1-  in  (46)  and  applying  Abel's  theorem. 


Corollary  4 

The  mean  L of  the  probability  density  {x^}  is  given  by 

1 -1  m-1 

(49)  L = -(m-a*)  { [a* (2)-m(m-l)+2m(m-a*) ) 11  - x^] 


+ x^ [ B*(2)-v(v-l)+2v(m-a*) ] } , 

where  a*(2)  and  B*(2),  0 < v < m-1,  are  the  second  factorial  moments  of  the 

densities  {a.}  and  (b  0 < v < m-1,  respectively. 

J vj  - - 


Proof 


By  differentiating  in  Formula  (26)  and  using  de  l'Hopital's  rule,  we 


obtain 


(50)  2(m-a*)L  = [a*(2)  - m(m-l)][l  - ^ xj  + ^ 6*(2)xv 

m-1  m-1  m-1 

£ v(v-l)x  + 2m  £ B*x  - 2a*  T.  vx  , 
v=0  v v=0  v v v=*0  v 

m-1 

and  upon  replacing  B*x^  by  the  expression,  obtained  from  (48)  , we  find 

the  expression  given  in  (49). 
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4.  Bailey's  Bulk  Queue 

One  of  the  earliest  papers  to  use  the  Rouche  method  of  solution  was 
N.  T.  J.  Bailey's  treatment  [1]  of  a bulk  queueing  model,  involving  a server 
who  becomes  available  at  the  epochs  of  a renewal  process  with  underlying 
distribution  H(’)-  Customers  arrive  according  to  a Poisson  process  of  rate  X. 
If  k customers  are  present  when  the  server  becomes  available,  a group  of  size 
min(k,m)  enters  service. 

It  is  well-known  that  the  successive  queue  lengths  immediately  prior  to 
the  beginnings  of  services  form  a Markov  chain  of  the  type  (1)  in  which 
b = a^,  for  0 < v < m-1,  and  j > 0.  The  quantities  a^  are  given  by 


,°°  -Xu  (Xu) J , tI/  s 

aj  = ;oe  ~jr~ dH(u)’ 


for  j > 0. 


We  shall  not  concern  ourselves  with  the  particular  form  (51),  but  it  is 
interesting  to  note  the  consequences  of  the  fact  that  the  first  m rows  are 
identical  and  given  by  {a^}. 

Theorem  5 

For  the  Markov  chain  in  Bailey's  model,  the  matrix  K(z)  defined  in 
Formula  (37),  has  m identical  rows  which  are  equal  to  the  first  row  of 
G(z) , defined  in  Formula  (13).  The  vector  jrt  defined  in  (30),  is  given  by 
the  first  row  of  the  matrix  G.  The  vector  _0,  defined  in  (31),  is  given  by 


1 = e, 


where  is  the  first  component  of  the  vector  £.  The  vector  C.  is  given  by 


L = hi  1 . 


and 
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(54) 


m-1 

E x = £ e = 
v=0  v 


The  mean  L is  given  by 

(55)  L = a*  + -|(m-a*)  1[a*(2)  + a*  + m2(p  1 - 1)  - E v2x  ]. 

2 1 v=0  v 

Proof 

' T 

It  is  trivial  to  verify  that  K(z)  has  identical  rows.  Comparing  the 
equation  (37)  and  the  equation  (14)  for  the  first  rows  only,  we  see  that 
these  equations  are  the  same.  This  shows  that  the  rows  of  K(z)  are  equal  to 
the  first  row  of  G(z).  The  statement  regarding  x is  now  obvious.  It 
follows  that  all  components  of  the  vector  6^  must  be  equal.  Comparing  the 
equations  (21)  and  (31)  for  the  first  components  of  and  0_,  we  obtain 
0^  = p^,  so  that  the  equation  (52)  follows.  Formula  (53)  follows  from 
Theorem  3 and  the  expression  (55)  for  L follows  by  substituting  (54)  in 
(49)  and  simplifying. 

A Markov  chain  of  the  type  (1)  occurs  also  in  the  study  of  the  stationary 

waiting  time  distribution  of  a discrete  version  of  the  GI/G/1  queue,  as 

investigated  by  G.  Ponstein  [14].  He  considers  service  and  interarrival 

distributions  of  lattice  type.  Choosing  the  step  of  the  lattice  as  unit  of 

time,  it  is  assumed  that  the  probability  density  of  S-T  (service  time 

minus  interarrival  time)  concentrates  on  the  integers  {-m, -m+l , . . . ,0, . . . } and 

that  P{S-T  = -m}  > 0.  Setting  P{S-T  = 1 ) = a . , , for  j > -m,  Ponstein  shows 

nri-j 

that  the  stationary  density  of  the  waiting  time  of  the  n-th  customer  is 
simply  related  to  the  stationary  density  of  a Markov  chain  of  the  type  which 
arises  in  Bailey's  model.  We  refer  to  [14]  for  the  details  of  an  algorithm, 
which  is  essentially  that  obtained  from  the  Rouche  method.  A substantial 
amount  of  work  is  devoted  to  the  problems  of  multiple  roots  and  in  some  cases, 


I 


see  Ponstein's  discussion  in  Section  6 of  [14],  numerical  difficulties  need 
to  be  overcome  by  tr ial-and-error . 

Our  preceding  discussion  provides  an  alternative  to  Ponstein's  algorithm. 
It  avoids  the  numerical  problems  of  his  method  and  appears  to  have  also  much 
smaller  computation  times. 


5.  Moran's  Dam  with  Infinite  Capacity 

A classical  model,  due  to  Moran  [5],  for  a dam  in  discrete  time  with 
discretized  content,  involves  a Markov  chain  of  the  type  (1)  in  which  the  first 
m rows  satisfy 


),_  = a , b = a,  .,  for  k > 1, 
jO  v=0  v jk  k+m-j 


and  0 < j < m-1.  In  this  model,  (a^)  is  the  probability  density  of  the 
number  of  units  of  water  added  to  the  dam  per  year  and  m is  the  maximum 
amount  of  water  released  at  the  end  of  each  year.  We  assume  that  the  capacity 
of  the  dam  is  infinite.  A detailed  discussion  of  this  model,  both  for  the 
finite  and  infinite  capacity  dam  is  given  in  Prabhu  [15],  Chapter  6.  In 
Moran's  model,  the  dam  content  is  described  immediately  following  releases  of 


water.  It  is  however  easy  to  see  that  if  we  consider  the  Bailey  queue, 
discussed  in  Section  4,  immediately  after  the  epochs  of  availability  of  the 
server,  we  obtain  precisely  the  embedded  Markov  chain  of  the  Moran  dam. 
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where  {x^}  is  the  density  obtained  in  Section  4.  The  expected  (stationary) 
dam  content  L',  following  releases  is  related  to  the  quantity  L of  Formula 


(55)  by 


- CO  - UI  -L  -V 

L'  = = L - m + (m-j)x^  = L - a*, 


since,  by  Formula  (48)  , we  have 


(59)  jIq  (m-j)Xj  = m - a*. 

We  see  that  Formula  (59)  also  gives  us  the  expected  deficiency  per  year  for 
the  stationary  dam. 

Let  q(mr+j-l) , for  r > 0,  1 < j < m,  be  the  expected  number  of 
periods  (years)  until  the  first  time  the  dam  is  deficient,  given  that  the 
dam  content  at  time  0+  is  mr+j-1,  then  we  have 


Theorem  6 

The  mean  q(mr+j-l)  is  given  by  the  j-th  component  of  the  vector 


[I  - Gr+1  + (r+l)G](I  - G + G)-1p  = y(r+1) 


Proof 


The  time  until  the  first  deficiency,  starting  with  a dam  content  of 
mr+j-1,  is  also  the  first  passage  time  from  the  state  (r+l,j)  to  the  level  0 
in  a Markov  chain  of  type  (1).  It  follows  that  the  vector  is  given  by 

,, , (r+1)  ,d  ~r+l.  ..  5 v 

(61.)  U = {-T-  G (z)}  e = G u, 

— dz  z=l  — v=0 


and  since 


ZQ  G ( I-G+G)  = I 


Gr+1  + (r+l)G, 


for  r > 0, 
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we  obtain  the  stated  result. 


Remarks 


a.  The  vectors  are  most  easily  computed  by  the  recurrence  relation 


(63) 


u(1)  = U,  p(r+1)  = £ + G p(r),  for  r > 1. 


b.  By  using  formulas  proved  in  [10] , it  is  possible  to  evaluate  also  the 
second,  and  in  principle,  higher  moments  of  this  waiting  time.  The  resulting 
expressions  are  however  quite  complicated  and  involve  the  matrix  M. 

c.  We  note  that  n(mr+j-l)  is  not  the  expected  time  till  the  first  emptiness, 
but  the  expected  number  of  years  until  the  dam  has  a positive  deficiency.  The 
formula  for  the  expected  time  until  the  first  emptiness  is  more  complicated, 
but  may  be  routinely  calculated. 

Corollary  5 

In  a dam  with  content  mr+j-1  at  time  0+,  the  expected  amount  x(mr+j-l) 
of  water  released  until  the  year  of  the  first  deficiency  is  given  by  the  j-th 
component  of  the  vector 


(64) 


(r+1)  ™ ,.  . . r+1 

my.  + i|1(i-l)  (G 


Proof 


The  first  term  is  obvious.  The  second  term  follows  by  recalling  that 


.r+1. 


(G  ) ^ is  the  conditional  probability  that  the  first  deficient  dam  content 
is  of  size  i-1,  or  equivalently  that  the  first  visit  to  level  0 occurs  in  the 
state  (0,i). 


V 
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6.  A Bulk  Service  Queue  with  a Threshold 

We  shall  now  consider  a variant  of  the  M/G/l  queue,  which  has  received 
wide  attention  (3,  6,  7].  This  model  has  a variety  of  concrete  applications, 
which  usually  involve  some  interesting  control  aspects.  We  shall  obtain  a 
number  of  explicit  and  computationally  tractable  results  not  available  before. 

Customers  arrive  at  a service  unit  according  to  a Poisson  process  of 
rate  X.  Services  occur  in  groups,  with  the  group  size  dependent  on  the  queue 
length  according  to  the  following  rule.  Let  there  be  i customers  waiting  at 
the  completion  of  a service.  If  0 < i < L,  the  server  remains  idle  until 
the  queue  length  reaches  i and  then  starts  serving  all  L customers.  If 
L < i < m,  a group  of  size  i enters  service  and  if  i > m,  a group  of  size  m 
is  served.  It  is  assumed  that  the  lengths  of  service  of  successive  groups 
are  conditionally  independent,  given  the  group  sizes  and  that  the  probability 
distribution  of  the  service  times  of  a group  of  size  j is  H^(*)  with  mean 
ctj , L < j < m.  For  notational  purposes,  we  denote  the  Laplace-Stielt jes 
transform  of  Hj(‘)  by  h^(s),  Re  s > 0. 

As  discussed  in  [7],  the  bivariate  sequence  of  the  queue  lengths 
following  departures  and  the  times  between  departures  defines  a Markov 
renewal  sequence,  whose  transition  probability  matrix  Q(’)  has  a Laplace- 
Stieltjes  transform  Q*(s),  which  is  of  the  form 


(65) 


b00(s) 

b01(s) 

b02(s) 

• * * 

Vs) 

bn  (s) 

b12(s) 

• • • 

Vl,0(,> 

b.-i,i<s) 

bm-l,2(s) 

• • • 

V’ 

a1(s) 

a2(s) 

• • • 

0 

a0(s) 

a^s) 

. . . 

0 

0 

a0(s) 

• • • 

• 

• 

• 

• 
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The  entries  b^(s),  0 < v < m-1,  j > 0,  are  given  by 


(66) 


bvj(s)  = (A^  f0  e"(A+S)u  dHL(u),  for  0 < v < L,  j > 0, 

for  L < v < m-1,  j > 0, 


and  the  entries  a^Cs)  are  given  by 


(67) 


/ , -(A+s)u  (Au)J  s 

a . (s)  = /_  e . dH  (u), 


for  j > 0. 


The  matrix  P = Q*(0+)  is  clearly  of  the  type  (1)  and 


(68) 


A*(z)  = h (A-Az) , 
m 


B^(z)  = hL(A-Az) , 


hv(A-Az), 


for  0 < v < L, 
for  L < v < m-1. 


a*  = A a 

m 


6*  = A aT  , 
v L 


= A a , 

v 


for  0 < v < L, 
for  L < v < m-1. 


Provided  A < m,  the  queue  is  stable  and  the  stationary  density  {x^}  of  the 
queue  length  following  departures  exists.  The  quantities  x^  can  be  computed 
by  using  the  procedure  discussed  in  Section  3. 


A.  The  Queue  Length  in  Continuous  Time 

In  this  section  we  investigate  the  stationary  density  of  the  queue 

length  at  an  arbitrary  point  of  time.  Henceforth  we  shall  only  consider  the 

stable  case  p*  = Am  <1. 

m 

We  shall  need  to  compute  the  fundamental  mean  E and  the  mean  recurrence 
times  Ej  of  the  states  j >0,  in  the  embedded  Markov  renewal  process  of  the 
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queue.  The  fundamental  mean  Is  the  inner  product  of  the  vector  x and  the 
vector  of  the  row  sum  means  of  the  transition  matrix  Q(.).  From  (65)  and 
(66)  we  obtain  that 


(69) 


E = 


IV 

A v^O 


i m-1 

(L-v)x  + — I 6*  x 
v A v=0  v v 


1 m-1 

+ T[1  - E.x  ]a*. 
A v=0  v 


Lemma  4 

The  fundamental  mean  E is  also  given  by 


(70) 


L-l 


m-1 


m-1 


E = — E x + — [1  - E x ] + -r  Z vx  , 
A v=0  v A v=0  v A v=L  v 


and  the  mean  recurrence  time  E.  of  the  state  j in  the  Markov  renewal  process 
Q(’)  is  given  by 


(71)  E.  = — , for  i > 0. 

IX.  J - 

J J 

Proof 

Formula  (48)  yields  that 


(72) 


m-1  m-1  m-1 

E 6*  x + a*[l  - Z x ] = m[l  - Ex] 
v=0  v v v=0  v v=0  v 


m-1 

+ E vx  , 
v=0  v 


and  substitution  in  (69)  leads  to  (70).  Formula  (71)  is  obtained  by  appli- 
cation of  Thm.  6.12,  p.  155  in  (Jinlar  [2]  or  Thm.  2.11,  p.  196  in  Hunter  [4]. 

We  shall  now  denote  the  density  of  the  queue  length  in  continuous  time 
by  2.  = That  density  is  related  to  the  density  x as  follows: 

Theorem  7 


The  density  is  given  by 


yj  XE  v=0  Xv’ 

1 L~1  f“>  -Xu  (Xu)  ^ L . u , n , 

yj  = ¥ vio  Vo  e -(fZTT  11  " HL(u)]du 


for  0 < j < L-l 


+ E v=L  Vo  e_XU  (-rr-_v>_!  U ' Vu)]du’  for  L < u < m-1 


1 r°°  “XU  (Xu)  ^ L r,  ,,  , s , , 

yj  ¥ v=o  Vo  e -JFl) T 11  »l  u du 


(j-L)! 

. 1 m"1  r°°  ~XU  (Xu)  d V , , ( t / \ l i 

+ — E,  x / e —7~. r-r  1 - H (u)Jdu 

E v=L  v 0 (j-v) ! v 


+ -jr  I x / e -■  V-yy  [1  - H (u)  ]du,  for  v > m. 
E v=m  v 0 ( j — v) ! m 


The  probability  generating  function 


Y(Z)  = jZ0  y.z  , 


for  I z | < 1, 


is  given  by 


1 k-±  \)  L 

— —r-j r { E X [z  - z h (X-Xz)] 

XE(l-z)  v=0  v L 


+ E x z [1  - h (X-Xz) ] 
v=L  v V 


m-1  v 

+ [X(z)  - v|Q  x^z  ][1  - hm(X-Xz)]}, 


The  means  L and  L*  of  x and  y_  are  related  by 


76)  2XEL*  = 2m(L  - yE  vx^)  + L(L-l)  J. 

m-1 

+ m[ (m-1)  - 2 (m-Xa  ) ] [1  - E.  x ] 
m v— u v 


+ 2XLa,  E.  x + 2X  E va  x , 
L v=0  v v=L  v v 
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m-1  2 (2)  L-l  2 m_l  (2) 

(77)  2 ( m-  A a ) (L  - E.  vx  ) = A a,  E_  x + A Er  av  'x 

' ' m V = 0 V * \ ) = ( 1 \>  \»=l  \» 


L v=0  v v=L  v v 


ro- 1 o ( ? A m—  1 

E v(v-l)x  + [A  oi  - m(m-l)  + 2m(m  - Aa  )][1  - E x], 
v=0  v m m v=0  v 


where  the  quantities  , L < v < m,  are  the  second  moments  of  the  probability 

distributions  H (•)• 


Proof 

The  formulas  (73)  were  proved  in  Neuts  [7]  by  a classical  argument  based 
on  the  Key  renewal  theorem.  The  probability  generating  function  Y(z)  is  ob- 
tained from  (73)  by  routine  but  lengthy  calculations  involving  a few  inter- 
changes of  summations  and  of  integrations.  Formula  (77)  is  obtained  by 
appropriate  substitutions  in  (49)  and  the  expression  (76)  is  obtained  by 
lengthy  differentiations  and  passage  to  the  limit  in  (75). 


Remark 

When  the  probability  distributions  H^O) , L < v < m,  are  of  phase  type 
[9],  the  integrals  which  appear  in  the  formulas  (73)  can  be  computed 
recursively  for  each  v.  No  numerical  integrations  are  then  required  and  the 
evaluation  of  the  quantities  y , j > 0,  from  the  x^,  j >0,  is  particularly 


Corollary  6 

In  the  stationary  version  of  the  queue  the  probability  Pq,  that  the 
server  is  idle  is  given  by 

L-l 

E (L-v)  x 

. . v=0  v 

0 " L=1  /T  . M m-1  . , m-1 

E-  (L-v)  x + m[l  - En  x ] + E_  vx 
V=0  v v=0  v v=0  v 
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